function g = grav_metric_eval(u, udt, udtdt, params)

  Nx = params.Nx;
  Ny = params.Ny;
  Nz = params.Nz;
  Nt = Nx.*Ny;

  g = zeros(Nx, Ny, Nz, 4, 4);


  T  = reshape(u(0.*Nt+1 : 1.*Nt), Nx, Ny);
  ux = reshape(u(1.*Nt+1 : 2.*Nt), Nx, Ny);
  uy = reshape(u(2.*Nt+1 : 3.*Nt), Nx, Ny);
  %ut = reshape(u(3*Nt+1 : 4*Nt), Nx, Ny);


  uxdt = reshape(udt(1.*Nt+1 : 2.*Nt), Nx, Ny);
  uydt = reshape(udt(2.*Nt+1 : 3.*Nt), Nx, Ny);
  %utdt = reshape(udt(3*Nt+1 : 4*Nt), Nx, Ny);

  uxdtdt = reshape(udtdt(1.*Nt+1 : 2.*Nt), Nx, Ny);
  uydtdt = reshape(udtdt(2.*Nt+1 : 3.*Nt), Nx, Ny);
  %utdtdt = reshape(udtdt(3*Nt+1 : 4*Nt), Nx, Ny);


  ut = sqrt(1+ux.^2+uy.^2);
  utdt = (uxdt.*ux + uydt.*uy)./ut;
  utdtdt = (1+ux.^2+uy.^2).^(-3./2).*(uxdt.^2.*(1+uy.^2)+(-2).*ux.*uxdt.*uy.*uydt+(1+ux.^2).*uydt.^2+(1+ux.^2+uy.^2).*(ux.*uxdtdt+uy.*uydtdt));

  utdx = dx(ut);
  utdy = dy(ut);

  uxdx = dx(ux);
  uxdy = dy(ux);

  uydx = dx(uy);
  uydy = dy(uy);


  utdtdx  = dx(utdt);
  utdtdy  = dy(utdt);
  utdxdx  = dx(utdx);
  utdxdy  = dy(utdx);
  utdydy  = dy(utdy);

  uxdtdx  = dx(uxdt);
  uxdtdy  = dy(uxdt);
  uxdxdx  = dx(uxdx);
  uxdxdy  = dy(uxdx);
  uxdydy  = dy(uxdy);

  uydtdx  = dx(uydt);
  uydtdy  = dy(uydt);
  uydxdx  = dx(uydx);
  uydxdy  = dy(uydx);
  uydydy  = dy(uydy);


  % g -- metric
  % gam -- gamma
  % gamd -- gamma derivatives
  global nh2s7 nk2s7 na2t2 na2t4 na2t5;
  if params.metricLinearFnInterpolation
    global interpgrid;
  end



  % *************************************************************
  % *************************************************************
  %   precompute x-y quantities before the r loop
  % *************************************************************
  % *************************************************************

  % all of the computations outside of the metric loops happen in terms of uU
  
  % *************************************************************
  % zeroth order
  % *************************************************************

  % *************************************************************
  % first order
  % *************************************************************

  % partial_lam u^lam
  delu = utdt+uxdx+uydy;

  % j_mu
  j1d = zeros(Nx, Ny, 3);
  j1d(:,:,1) = -(ut.*utdt+ux.*utdx+uy.*utdy);
  j1d(:,:,2) =  (ut.*uxdt+ux.*uxdx+uy.*uxdy);
  j1d(:,:,3) =  (ut.*uydt+ux.*uydx+uy.*uydy);


  % sigma_\mu_\nu

  a1dd = zeros(Nx, Ny, 3, 3);
  a1dd(:,:,1,1) = (1/2).*((-1).*((-1)+ut.^2).^2.*utdt+uxdx+ut.*(utdx.*ux+(-1).*ux.*uxdt+utdy.*uy+(-1).*uy.*uydt)+ut.^3.*((-1).*utdx.*ux+ux.*uxdt+(-1).*utdy.*uy+uy.*uydt)+ux.*(ux.*uxdx+uy.*(uxdy+uydx))+uydy+uy.^2.*uydy+ut.^2.*(((-1)+ux.^2).*uxdx+ux.*uy.*(uxdy+uydx)+((-1)+uy.^2).*uydy));
  a1dd(:,:,1,2) = (1/2).*(ut.^3.*utdt.*ux+((-1)+ut.^2).*utdx.*(1+ux.^2)+uxdt+ux.*(ux.*uxdt+uy.*((-1).*utdy+uydt))+(-1).*ut.^2.*((1+ux.^2).*uxdt+ux.*uy.*((-1).*utdy+uydt))+(-1).*ut.*(utdt.*ux+(1+ux.^2).*(ux.*uxdx+uy.*(uxdy+uydx))+ux.*((-1)+uy.^2).*uydy));
  a1dd(:,:,1,3) = (1/2).*(ut.^3.*utdt.*uy+((-1)+ut.^2).*utdy.*(1+uy.^2)+uydt+uy.*((-1).*utdx.*ux+ux.*uxdt+uy.*uydt)+ut.^2.*(ux.*(utdx+(-1).*uxdt).*uy+(-1).*(1+uy.^2).*uydt)+(-1).*ut.*(ux.^2.*uxdx.*uy+ux.*(1+uy.^2).*(uxdy+uydx)+uy.*(utdt+(-1).*uxdx+uydy+uy.^2.*uydy)));
  a1dd(:,:,2,2) = (1/2).*(utdt.*((-1)+ut.^2+(-1).*(1+ut.^2).*ux.^2)+ut.*((-1).*(ux+ux.^3).*(utdx+(-1).*uxdt)+utdy.*uy+(-1).*utdy.*ux.^2.*uy+((-1)+ux.^2).*uy.*uydt)+(1+ux.^2).*((1+ux.^2).*uxdx+ux.*uy.*(uxdy+uydx))+((-1)+(-1).*uy.^2+ux.^2.*((-1)+uy.^2)).*uydy);
  a1dd(:,:,2,3) = (1/2).*(uxdy+(-1).*ut.^2.*utdt.*ux.*uy+ux.^3.*uxdx.*uy+ut.*((-1).*(1+ux.^2).*(utdx+(-1).*uxdt).*uy+(-1).*utdy.*ux.*(1+uy.^2)+ux.*(1+uy.^2).*uydt)+uydx+uy.^2.*(uxdy+uydx)+ux.^2.*(1+uy.^2).*(uxdy+uydx)+ux.*uy.*((-1).*utdt+uxdx+uydy+uy.^2.*uydy));
  a1dd(:,:,3,3) = (1/2).*((-1).*uxdx+(-1).*ux.^2.*uxdx+ux.*uxdy.*uy+(-1).*uxdx.*uy.^2+ux.^2.*uxdx.*uy.^2+ux.*uxdy.*uy.^3+utdt.*((-1)+ut.^2+(-1).*(1+ut.^2).*uy.^2)+(-1).*ut.*(ux.*((-1).*utdx+uxdt)+utdy.*uy+ux.*(utdx+(-1).*uxdt).*uy.^2+utdy.*uy.^3)+ut.*(uy+uy.^3).*uydt+ux.*uy.*uydx+ux.*uy.^3.*uydx+(1+uy.^2).^2.*uydy);


  % *****************************************************
  %  second order x^\mu dependent quantities
  % *****************************************************

  % TODO all of the scalar, vector and tensor terms are in nathan's mathematica file, but somebody needs to check them (rederive?) should be a 30min job
  % TODO scalar terms
  s3 = utdt.^2+(-1).*uxdt.^2+(-1).*uydt.^2+ut.*(utdtdt+uxdtdx+uydtdy)+ux.^2.*(utdt.^2+utdx.^2+(-1).*uxdt.^2+(-1).*uxdx.^2+(-1).*uydt.^2+(-1).*uydx.^2)+ux.*(utdtdx+uxdxdx+ut.*(2.*utdt.*utdx+(-2).*uxdt.*uxdx+(-2).*uydt.*uydx)+uydxdy)+uy.^2.*(utdt.^2+utdy.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(-1).*uydt.^2+(-1).*uydy.^2)+uy.*(utdtdy+uxdxdy+ut.*(2.*utdt.*utdy+(-2).*uxdt.*uxdy+(-2).*uydt.*uydy)+ux.*(2.*utdx.*utdy+(-2).*uxdx.*uxdy+(-2).*uydx.*uydy)+uydydy);
  s4 = (-1).*utdt.^2+uxdt.^2+uydt.^2+2.*ut.*ux.*((-1).*utdt.*utdx+uxdt.*uxdx+uydt.*uydx)+ux.^2.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+uy.^2.*((-1).*utdt.^2+(-1).*utdy.^2+uxdt.^2+uxdy.^2+uydt.^2+uydy.^2)+uy.*(2.*ut.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+2.*ux.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy));
  s5 = (utdt+uxdx+uydy).^2;
  s6 = ((-1).*(utdx+uxdt).*uy+ux.*(utdy+uydt)+ut.*((-1).*uxdy+uydx)).^2;
  s7 = (1/2).*(uxdx.^2+(-2).*ut.*ux.^3.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+ux.^4.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2)+uxdy.^2+2.*uxdy.*uydx+uydx.^2+(-2).*ut.*ux.^2.*uy.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+2.*ux.^3.*uy.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx))+(-2).*uxdx.*uydy+uydy.^2+(-2).*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+2.*ux.*uy.*(uxdy+uydx).*((-2).*utdt+uxdx+uydy)+(-2).*ut.*ux.*uy.^2.*((utdy+(-1).*uydt).*(uxdy+uydx)+(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+ux.^2.*uy.^2.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy))+2.*ux.*uy.^3.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+2.*ut.*ux.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(utdx+(-1).*uxdt).*((-1).*uxdx+uydy))+ux.^2.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+(-2).*((-1).*utdt+uxdx).*((-1).*uxdx+uydy))+(-2).*ut.*uy.*((utdx+(-1).*uxdt).*(uxdy+uydx)+(utdy+(-1).*uydt).*((-1).*uxdx+uydy))+uy.^2.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uydy).*((-1).*uxdx+uydy)));


  v3d = zeros(Nx, Ny, 3);
  v4d = zeros(Nx, Ny, 3);
  v5d = zeros(Nx, Ny, 3);
  v6d = zeros(Nx, Ny, 3);
  v7d = zeros(Nx, Ny, 3);

  % TODO check the vector expressions. 
  v3d(:,:,1) = ux.*((-1).*utdt.*utdx.*(1+2.*ux.^2)+(1+2.*ux.^2).*uxdt.*uxdx+(-1).*ux.*(utdtdt+uxdtdx+uydtdy)+uydt.*uydx+2.*ux.^2.*uydt.*uydx)+(-1).*uy.^2.*(utdtdt+uxdtdx+uydtdy+ux.*(2.*utdt.*utdx+(-2).*(uxdt.*uxdx+uydt.*uydx)))+(-1).*(1+2.*ux.^2).*uy.*(utdt.*utdy+(-1).*uxdt.*uxdy+(-1).*uydt.*uydy)+2.*uy.^3.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+ut.*(ux.^2.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+(-1).*ux.*(utdtdx+uxdxdx+uydxdy)+2.*ux.*uy.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy)+uy.^2.*((-1).*utdt.^2+(-1).*utdy.^2+uxdt.^2+uxdy.^2+uydt.^2+uydy.^2)+(-1).*uy.*(utdtdy+uxdxdy+uydydy));
  v3d(:,:,2) = utdtdx+uxdxdx+ut.*(utdt.*utdx+(-1).*uxdt.*uxdx+(-1).*uydt.*uydx)+ux.^3.*(utdt.^2+utdx.^2+(-1).*uxdt.^2+(-1).*uxdx.^2+(-1).*uydt.^2+(-1).*uydx.^2)+uydxdy+uy.*(utdx.*utdy+(-1).*uxdx.*uxdy+(-1).*uydx.*uydy)+ux.^2.*(utdtdx+uxdxdx+ut.*(2.*utdt.*utdx+(-2).*uxdt.*uxdx+(-2).*uydt.*uydx)+uydxdy+uy.*(2.*utdx.*utdy+(-2).*uxdx.*uxdy+(-2).*uydx.*uydy))+ux.*(utdt.^2+utdx.^2+(-1).*uxdt.^2+(-1).*uxdx.^2+(-1).*uydt.^2+ut.*(utdtdt+uxdtdx+uydtdy)+(-1).*uydx.^2+uy.^2.*(utdt.^2+utdy.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(-1).*uydt.^2+(-1).*uydy.^2)+uy.*(utdtdy+uxdxdy+ut.*(2.*utdt.*utdy+(-2).*uxdt.*uxdy+(-2).*uydt.*uydy)+uydydy));
  v3d(:,:,3) = utdtdy+uxdxdy+ut.*(utdt.*utdy+(-1).*uxdt.*uxdy+(-1).*uydt.*uydy)+ux.*(utdx.*utdy+(-1).*uxdx.*uxdy+(-1).*uydx.*uydy)+uy.^3.*(utdt.^2+utdy.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(-1).*uydt.^2+(-1).*uydy.^2)+uy.*(utdt.^2+utdy.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(-1).*uydt.^2+ut.*(utdtdt+uxdtdx+uydtdy)+ux.^2.*(utdt.^2+utdx.^2+(-1).*uxdt.^2+(-1).*uxdx.^2+(-1).*uydt.^2+(-1).*uydx.^2)+ux.*(utdtdx+uxdxdx+ut.*(2.*utdt.*utdx+(-2).*uxdt.*uxdx+(-2).*uydt.*uydx)+uydxdy)+(-1).*uydy.^2)+uydydy+uy.^2.*(utdtdy+uxdxdy+ut.*(2.*utdt.*utdy+(-2).*uxdt.*uxdy+(-2).*uydt.*uydy)+ux.*(2.*utdx.*utdy+(-2).*uxdx.*uxdy+(-2).*uydx.*uydy)+uydydy);


  v4d(:,:,1) = ux.^4.*(utdtdt+utdxdx+(-2).*uxdtdx)+ux.^3.*((-1).*ut.*((-2).*utdtdx+uxdtdt+uxdxdx)+2.*uy.*(utdxdy+(-1).*uxdtdy+(-1).*uydtdx))+ux.^2.*(utdxdx+utdydy+(-2).*uxdtdx+uy.^2.*(utdxdx+utdydy+(-2).*((-1).*utdtdt+uxdtdx+uydtdy))+(-1).*ut.*uy.*((-2).*utdtdy+2.*uxdxdy+uydtdt+uydxdx))+ux.*((-1).*ut.*(uxdxdx+uxdydy)+2.*uy.^3.*(utdxdy+(-1).*uxdtdy+(-1).*uydtdx)+(-2).*uy.*(uxdtdy+uydtdx)+(-1).*ut.*uy.^2.*((-2).*utdtdx+uxdtdt+uxdydy+2.*uydxdy))+uy.*(uy.^3.*(utdtdt+utdydy+(-2).*uydtdy)+uy.*(utdxdx+utdydy+(-2).*uydtdy)+(-1).*ut.*uy.^2.*((-2).*utdtdy+uydtdt+uydydy)+(-1).*ut.*(uydxdx+uydydy));
  v4d(:,:,2) = uxdxdx+ux.^4.*((-2).*utdtdx+uxdtdt+uxdxdx)+uxdydy+(uxdtdt+uxdydy).*uy.^2+ut.*((-1).*ux.^3.*(utdtdt+utdxdx+(-2).*uxdtdx)+2.*uxdtdy.*uy+2.*ux.^2.*uy.*((-1).*utdxdy+uxdtdy+uydtdx)+ux.*((-1).*utdxdx+2.*uxdtdx+(-1).*utdydy.*(1+uy.^2)+uy.^2.*((-1).*utdtdt+2.*uydtdy)))+ux.^3.*uy.*((-2).*utdtdy+2.*uxdxdy+uydtdt+uydxdx)+ux.^2.*((-2).*utdtdx+uxdtdt+2.*uxdxdx+uxdydy+uy.^2.*((-2).*utdtdx+uxdtdt+uxdydy+2.*uydxdy))+ux.*uy.*(2.*uxdxdy+(-2).*utdtdy.*(1+uy.^2)+uy.^2.*uydtdt+uydxdx+(1+uy.^2).*uydydy);
  v4d(:,:,3) = ux.^3.*((-2).*utdtdx+uxdtdt+uxdxdx).*uy+ut.*(2.*ux.*uydtdx+2.*ux.*uy.^2.*((-1).*utdxdy+uxdtdy+uydtdx)+(-1).*uy.^3.*(utdtdt+utdydy+(-2).*uydtdy)+uy.*((-1).*utdxdx+(-1).*utdydy+(-1).*ux.^2.*(utdtdt+utdxdx+(-2).*uxdtdx)+2.*uydtdy))+ux.^2.*(uydtdt+uydxdx+uy.^2.*((-2).*utdtdy+2.*uxdxdy+uydtdt+uydxdx))+ux.*uy.*(uxdxdx+uxdtdt.*uy.^2+(-2).*utdtdx.*(1+uy.^2)+uxdydy.*(1+uy.^2)+2.*(1+uy.^2).*uydxdy)+(1+uy.^2).*(uydxdx+uydydy+uy.^2.*((-2).*utdtdy+uydtdt+uydydy));

  v5d(:,:,1) = (-1).*(utdt+uxdx+uydy).*((-1).*ut.*((-1)+ut.^2).*utdt+ux.*((-1).*((-1)+ut.^2).*utdx+ut.^2.*uxdt)+ut.*ux.^2.*uxdx+uy.*((-1).*((-1)+ut.^2).*utdy+ut.*(ut.*uydt+ux.*(uxdy+uydx)))+ut.*uy.^2.*uydy);
  v5d(:,:,2) = (utdt+uxdx+uydy).*(ux.*(1+ux.^2).*((-1).*utdt+uxdx)+ut.*(uxdt+ux.^2.*((-1).*utdx+uxdt)+ux.*uy.*((-1).*utdy+uydt))+uy.*(uxdy+ux.^2.*(uxdy+uydx))+ux.*uy.^2.*((-1).*utdt+uydy));
  v5d(:,:,3) = (utdt+uxdx+uydy).*(ux.^2.*((-1).*utdt+uxdx).*uy+ut.*(ux.*((-1).*utdx+uxdt).*uy+uydt+uy.^2.*((-1).*utdy+uydt))+ux.*(uydx+uy.^2.*(uxdy+uydx))+uy.*(1+uy.^2).*((-1).*utdt+uydy));

  v6d(:,:,1) = (1+2.*ux.^2+2.*uy.^2).*(ux.*(utdt.*utdx+(-1).*uxdt.*uxdx+(-1).*uydt.*uydx)+uy.*(utdt.*utdy+(-1).*uxdt.*uxdy+(-1).*uydt.*uydy))+ut.*(ux.^2.*(utdt.^2+utdx.^2+(-1).*uxdt.^2+(-1).*uxdx.^2+(-1).*uydt.^2+(-1).*uydx.^2)+ux.*uy.*(2.*utdx.*utdy+(-2).*uxdx.*uxdy+(-2).*uydx.*uydy)+uy.^2.*(utdt.^2+utdy.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(-1).*uydt.^2+(-1).*uydy.^2));
  v6d(:,:,2) = ut.*((-1).*utdt.*utdx+uxdt.*uxdx+uydt.*uydx)+ux.^3.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+uy.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy)+ux.^2.*(2.*ut.*((-1).*utdt.*utdx+uxdt.*uxdx+uydt.*uydx)+2.*uy.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy))+ux.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2+2.*ut.*uy.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+uy.^2.*((-1).*utdt.^2+(-1).*utdy.^2+uxdt.^2+uxdy.^2+uydt.^2+uydy.^2));
  v6d(:,:,3) = ut.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+ux.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy)+uy.^3.*((-1).*utdt.^2+(-1).*utdy.^2+uxdt.^2+uxdy.^2+uydt.^2+uydy.^2)+uy.*((-1).*utdt.^2+(-1).*utdy.^2+uxdt.^2+uxdy.^2+uydt.^2+2.*ut.*ux.*((-1).*utdt.*utdx+uxdt.*uxdx+uydt.*uydx)+ux.^2.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+uydy.^2)+uy.^2.*(2.*ut.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+2.*ux.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy));

  v7d(:,:,1) = ut.*((-1)+ut.^2).*(utdt.^2+utdx.*uxdt+utdy.*uydt)+ux.*(((-1)+ut.^2).*utdx.*(utdt+uxdx)+(-1).*ut.^2.*(utdt.*uxdt+uxdt.*uxdx+uxdy.*uydt)+((-1)+ut.^2).*utdy.*uydx)+(-1).*ut.*ux.^2.*(utdx.*uxdt+uxdx.^2+uxdy.*uydx)+(-1).*ut.*uy.^2.*(utdy.*uydt+uxdy.*uydx+uydy.^2)+uy.*((-1).*utdx.*uxdy+utdy.*(((-1)+ut.^2).*utdt+(-1).*ut.*ux.*uxdt+((-1)+ut.^2).*uydy)+ut.^2.*(utdx.*uxdy+(-1).*uxdt.*uydx+(-1).*uydt.*(utdt+uydy))+(-1).*ut.*ux.*(utdx.*uydt+uxdx.*uydx+uydx.*uydy+uxdy.*(uxdx+uydy)));
  v7d(:,:,2) = (-1).*ut.^2.*ux.*(utdt.^2+utdx.*uxdt+utdy.*uydt)+ux.*(1+ux.^2).*(utdx.*uxdt+uxdx.^2+uxdy.*uydx)+ux.*uy.^2.*(utdy.*uydt+uxdy.*uydx+uydy.^2)+uy.*(utdy.*uxdt+(1+ux.^2).*uxdy.*(uxdx+uydy)+ux.^2.*(utdy.*uxdt+utdx.*uydt+uxdx.*uydx+uydx.*uydy))+ut.*(utdt.*uxdt+uxdt.*uxdx+uxdy.*uydt+ux.^2.*(utdt.*uxdt+uxdt.*uxdx+(-1).*utdx.*(utdt+uxdx)+uxdy.*uydt+(-1).*utdy.*uydx)+ux.*uy.*((-1).*utdx.*uxdy+utdt.*uydt+uxdt.*uydx+uydt.*uydy+(-1).*utdy.*(utdt+uydy)));
  v7d(:,:,3) = (-1).*ut.^2.*uy.*(utdt.^2+utdx.*uxdt+utdy.*uydt)+ux.^2.*uy.*(utdx.*uxdt+uxdx.^2+uxdy.*uydx)+uy.*(1+uy.^2).*(utdy.*uydt+uxdy.*uydx+uydy.^2)+ut.*(uxdt.*uydx+ux.*uy.*(utdt.*uxdt+uxdt.*uxdx+(-1).*utdx.*(utdt+uxdx)+uxdy.*uydt+(-1).*utdy.*uydx)+uydt.*(utdt+uydy)+uy.^2.*((-1).*utdx.*uxdy+utdt.*uydt+uxdt.*uydx+uydt.*uydy+(-1).*utdy.*(utdt+uydy)))+ux.*(utdx.*uydt+uxdx.*uydx+uydx.*uydy+uy.^2.*(utdy.*uxdt+utdx.*uydt+uxdx.*uydx+uydx.*uydy+uxdy.*(uxdx+uydy)));


  %*********************************************
  % tensors


  t2dd = zeros(Nx, Ny, 3, 3);
  t3dd = zeros(Nx, Ny, 3, 3);
  t4dd = zeros(Nx, Ny, 3, 3);
  t5dd = zeros(Nx, Ny, 3, 3);

  t2dd(:,:,1,1) = (1/2).*(ux.^5.*((-2).*utdtdx+uxdtdt+uxdxdx)+ut.*((-1).*ux.^4.*(utdtdt+utdxdx+(-2).*uxdtdx)+2.*ux.^3.*uy.*((-1).*utdxdy+uxdtdy+uydtdx)+2.*ux.*uy.*(uxdtdy+uydtdx+uy.^2.*((-1).*utdxdy+uxdtdy+uydtdx))+uy.^2.*((-1).*uxdtdx+(-1).*uy.^2.*(utdtdt+utdydy+(-2).*uydtdy)+uydtdy)+ux.^2.*(uxdtdx+(-1).*uydtdy+(-1).*uy.^2.*(utdxdx+utdydy+(-2).*((-1).*utdtdt+uxdtdx+uydtdy))))+ux.^4.*uy.*((-2).*utdtdy+2.*uxdxdy+uydtdt+uydxdx)+ux.*uy.^2.*((-1).*uxdxdx+uxdtdt.*(1+uy.^2)+uxdydy.*(2+uy.^2)+(-1).*utdtdx.*(1+2.*uy.^2)+(3+2.*uy.^2).*uydxdy)+ux.^3.*((-1).*utdtdx+uxdtdt+uxdxdx+(-1).*uydxdy+uy.^2.*((-4).*utdtdx+2.*uxdtdt+uxdxdx+uxdydy+2.*uydxdy))+ux.^2.*uy.*((-1).*utdtdy+uxdxdy.*(3+2.*uy.^2)+uydtdt+2.*uydxdx+uy.^2.*((-4).*utdtdy+2.*uydtdt+uydxdx)+((-1)+uy.^2).*uydydy)+uy.^3.*((-1).*uxdxdy+(-1).*utdtdy.*(1+2.*uy.^2)+(1+uy.^2).*uydtdt+(1+uy.^2).*uydydy));
  t2dd(:,:,1,2) = -(1/2).*((-1).*ux.^5.*(utdtdt+utdxdx+(-2).*uxdtdx)+ux.^4.*(ut.*((-2).*utdtdx+uxdtdt+uxdxdx)+2.*uy.*((-1).*utdxdy+uxdtdy+uydtdx))+ux.^3.*((-1).*utdtdt+(-1).*utdxdx+3.*uxdtdx+(-1).*uydtdy+(-1).*uy.^2.*(utdxdx+utdydy+(-2).*((-1).*utdtdt+uxdtdx+uydtdy))+ut.*uy.*((-2).*utdtdy+2.*uxdxdy+uydtdt+uydxdx))+uy.*(uxdtdy+uydtdx+uy.^2.*((-1).*utdxdy+2.*uxdtdy+uydtdx)+ut.*uy.*((-1).*utdtdx+uxdtdt+uxdydy+uydxdy))+ux.^2.*(uy.*((-1).*utdxdy.*(1+2.*uy.^2)+uxdtdy.*(3+2.*uy.^2)+2.*(1+uy.^2).*uydtdx)+ut.*((-1).*utdtdx+uxdtdt+uxdxdx+(-1).*uydxdy+uy.^2.*((-2).*utdtdx+uxdtdt+uxdydy+2.*uydxdy)))+ux.*(uxdtdx+(-1).*(utdtdt+utdxdx+(-2).*uxdtdx).*uy.^2+(-1).*uy.^4.*(utdtdt+utdydy+(-2).*uydtdy)+(-1).*uydtdy+ut.*uy.*(2.*uxdxdy+uydxdx+(-1).*uydydy)+ut.*uy.^3.*((-2).*utdtdy+uydtdt+uydydy)));
  t2dd(:,:,1,3) = -(1/2).*((-1).*ux.^4.*(utdtdt+utdxdx+(-2).*uxdtdx).*uy+ux.^3.*(ut.*((-2).*utdtdx+uxdtdt+uxdxdx).*uy+(-1).*utdxdy.*(1+2.*uy.^2)+uxdtdy.*(1+2.*uy.^2)+2.*(1+uy.^2).*uydtdx)+ux.^2.*((-1).*uy.*(utdtdt+utdydy+(-2).*uydtdy)+(-1).*uy.^3.*(utdxdx+utdydy+(-2).*((-1).*utdtdt+uxdtdx+uydtdy))+ut.*((-1).*utdtdy+uxdxdy+uydtdt+uydxdx)+ut.*uy.^2.*((-2).*utdtdy+2.*uxdxdy+uydtdt+uydxdx))+ux.*(uxdtdy+uydtdx+uy.*(uy.*(2.*uxdtdy.*(1+uy.^2)+(-1).*utdxdy.*(1+2.*uy.^2)+(3+2.*uy.^2).*uydtdx)+ut.*((-1).*uxdxdx+((-2).*utdtdx+uxdtdt).*uy.^2+uxdydy.*(1+uy.^2)+2.*(1+uy.^2).*uydxdy)))+uy.*((-1).*uxdtdx+(-1).*uy.^2.*(utdtdt+utdydy+uxdtdx+(-3).*uydtdy)+(-1).*uy.^4.*(utdtdt+utdydy+(-2).*uydtdy)+uydtdy+ut.*uy.^3.*((-2).*utdtdy+uydtdt+uydydy)+ut.*uy.*((-1).*utdtdy+(-1).*uxdxdy+uydtdt+uydydy)));
  t2dd(:,:,2,2) = (1/2).*(ux.^5.*((-2).*utdtdx+uxdtdt+uxdxdx)+ut.*(2.*ux.*uy.*(uxdtdy+ux.^2.*((-1).*utdxdy+uxdtdy+uydtdx))+(-1).*((-1)+ux.^2).*uy.^2.*(utdtdt+utdydy+(-2).*uydtdy)+(-1).*(1+ux.^2).*(ux.^2.*(utdtdt+utdxdx+(-2).*uxdtdx)+(-1).*uxdtdx+uydtdy))+ux.^4.*uy.*((-2).*utdtdy+2.*uxdxdy+uydtdt+uydxdx)+ux.*((-1).*utdtdx+uxdtdt+uxdxdx+(uxdtdt+uxdydy).*uy.^2+(-1).*uydxdy)+ux.^3.*((-3).*utdtdx+2.*uxdtdt+2.*uxdxdx+(-1).*uydxdy+uy.^2.*((-2).*utdtdx+uxdtdt+uxdydy+2.*uydxdy))+ux.^2.*uy.*(3.*uxdxdy+(-1).*utdtdy.*(1+2.*uy.^2)+uy.^2.*uydtdt+uydxdx+((-1)+uy.^2).*uydydy)+uy.*(uxdxdy+utdtdy.*(1+2.*uy.^2)+(-1).*(1+uy.^2).*uydtdt+(-1).*(1+uy.^2).*uydydy));
  t2dd(:,:,2,3) = (1/2).*(ut.*(uxdtdy+(-1).*ux.^3.*(utdtdt+utdxdx+(-2).*uxdtdx).*uy+uydtdx+uy.^2.*((-1).*utdxdy+2.*uxdtdy+uydtdx)+ux.^2.*((-1).*(utdxdy+(-1).*uxdtdy).*(1+2.*uy.^2)+2.*(1+uy.^2).*uydtdx)+ux.*uy.*((-1).*utdxdx+(-1).*utdydy.*(1+uy.^2)+(-1).*utdtdt.*(2+uy.^2)+2.*(uxdtdx+(1+uy.^2).*uydtdy)))+ux.*(1+ux.^2).*((-1).*utdtdy+uxdxdy+uydtdt+uydxdx)+uy.^3.*((1+ux.^2).*uxdtdt+(1+ux.^2).*uxdydy+(1+2.*ux.^2).*((-1).*utdtdx+uydxdy))+uy.*((-1).*utdtdx+uxdtdt+(1+ux.^2).*uxdydy+uydxdy+ux.^2.*((-1).*(2+ux.^2).*(2.*utdtdx+(-1).*uxdtdt)+(1+ux.^2).*uxdxdx+2.*uydxdy))+ux.*uy.^4.*((-2).*utdtdy+uydtdt+uydydy)+ux.*uy.^2.*(2.*(1+ux.^2).*uxdxdy+(-1).*(2+ux.^2).*(2.*utdtdy+(-1).*uydtdt)+(1+ux.^2).*uydxdx+uydydy));
  t2dd(:,:,3,3) = (1/2).*(ux.^3.*((-2).*utdtdx+uxdtdt+uxdxdx).*((-1)+uy.^2)+ut.*(ux.^2.*(utdtdt+utdxdx+(-2).*uxdtdx)+(-1).*uxdtdx+2.*ux.*uy.*uydtdx+2.*ux.*uy.^3.*((-1).*utdxdy+uxdtdy+uydtdx)+(-1).*uy.^4.*(utdtdt+utdydy+(-2).*uydtdy)+uydtdy+uy.^2.*((-1).*utdtdt+(-1).*utdydy+(-1).*ux.^2.*(utdtdt+utdxdx+(-2).*uxdtdx)+(-1).*uxdtdx+3.*uydtdy))+ux.^2.*uy.*(uydtdt+uydxdx+uy.^2.*((-2).*utdtdy+2.*uxdxdy+uydtdt+uydxdx))+ux.*(1+uy.^2).*(utdtdx+(-1).*uxdtdt+(-1).*uxdxdx+uydxdy+uy.^2.*((-2).*utdtdx+uxdtdt+uxdydy+2.*uydxdy))+uy.*(1+uy.^2).*((-1).*uxdxdy+(-1).*utdtdy.*(1+2.*uy.^2)+(1+uy.^2).*uydtdt+(1+uy.^2).*uydydy));
  
  t3dd(:,:,1,1) = (1/2).*(ux.^6.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2)+ux.^5.*((-2).*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+2.*uy.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+uy.^2.*((-1).*uxdt.^2+uydt.^2+(-2).*ut.*uy.*(uxdt.*uxdy+uydt.*(utdt+(-1).*uydy))+(-2).*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+uy.^2.*((-1).*uxdt.^2+(-1).*uxdy.^2+2.*uydt.*((-1).*utdy+uydt)+((-1).*utdt+uydy).^2))+ux.^3.*(2.*ut.*(uxdt.*((-1).*utdt+uxdx)+(-1).*uydt.*uydx)+(-2).*uy.*(utdy.*uxdt+(utdt+(-1).*uxdx).*uxdy+(utdx+(-3).*uxdt).*uydt+uydx.*(utdt+(-2).*uxdx+uydy))+(-2).*ut.*uy.^2.*((utdy+(-1).*uydt).*(uxdy+uydx)+(utdx+(-1).*uxdt).*((-2).*utdt+uxdx+uydy))+2.*uy.^3.*(2.*(utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-2).*utdt+uxdx+uydy)))+ux.^2.*(uxdt.^2+(-1).*uydt.^2+2.*ut.*uy.*(uxdt.*(uxdy+2.*uydx)+(-1).*uydt.*(utdt+(-2).*uxdx+uydy))+uy.^2.*(2.*utdt.^2+(-2).*utdx.*uxdt+uxdt.^2+(-1).*uxdx.^2+uxdy.^2+(-2).*utdy.*uydt+uydt.^2+4.*uxdy.*uydx+uydx.^2+4.*uxdx.*uydy+(-1).*uydy.^2+(-2).*utdt.*(uxdx+uydy))+2.*ut.*uy.^3.*((-1).*(utdx+(-1).*uxdt).*(uxdy+uydx)+(-1).*(utdy+(-1).*uydt).*((-2).*utdt+uxdx+uydy))+uy.^4.*((utdx+(-1).*uxdt).^2+2.*(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+((-1).*utdt+uydy).*((-3).*utdt+2.*uxdx+uydy)))+ux.^4.*(2.*uxdt.*((-1).*utdx+uxdt)+((-1).*utdt+uxdx).^2+(-1).*uydt.^2+(-1).*uydx.^2+(-2).*ut.*uy.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+uy.^2.*(2.*(utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+((-1).*utdt+uxdx).*((-3).*utdt+uxdx+2.*uydy)))+ux.*uy.*(4.*uxdt.*uydt+(-2).*uy.^2.*(utdy.*uxdt+(utdx+(-3).*uxdt).*uydt+uxdy.*(utdt+uxdx+(-2).*uydy)+uydx.*(utdt+(-1).*uydy))+2.*uy.^4.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+2.*ut.*uy.*(uydt.*(2.*uxdy+uydx)+uxdt.*((-1).*utdt+(-1).*uxdx+2.*uydy)+uy.^2.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+uydy)))));
  t3dd(:,:,1,2) = -(1/2).*(ut.*(ux.^5.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2)+2.*uxdt.*(uy+uy.^3).*uydt+2.*ux.^4.*uy.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx))+2.*uy.^3.*((-1).*utdy.*uxdt+uxdy.*((-1).*utdt+uydy))+ux.^3.*(2.*uxdt.*((-1).*utdx+uxdt)+((-1).*utdt+uxdx).^2+(-1).*uydt.^2+(-1).*uydx.^2+uy.^2.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy)))+2.*ux.^2.*uy.*(((-1).*utdt+uxdx).*uxdy+(-1).*uxdt.*(utdy+(-2).*uydt)+uydx.*(uxdx+(-1).*uydy)+uy.^2.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy)))+ux.*(uxdt.^2+(-1).*uydt.^2+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+uy.^2.*((-2).*utdx.*uxdt+uxdt.^2+uxdy.^2+2.*uxdy.*uydx+(-1).*((-1).*utdt+uydy).*(utdt+(-2).*uxdx+uydy))))+2.*((-1).*ux.^6.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+ux.^5.*uy.*((-1).*((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(-1).*(utdx+(-1).*uxdt).*(uxdy+uydx))+ux.*uy.*(uxdt.*(uxdy+uydx)+uy.^2.*((-1).*utdx.*uxdy+(-1).*((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+uxdt.*(uxdy+uydx))+uydt.*(uxdx+(-1).*uydy)+(-1).*uy.^4.*(utdy+(-1).*uydt).*((-1).*utdt+uydy))+uy.^2.*(uxdy.*uydt+uxdt.*((-1).*utdt+uydy)+uy.^2.*(uxdy.*((-1).*utdy+uydt)+uxdt.*((-1).*utdt+uydy)))+ux.^4.*((-1).*(utdx+(-2).*uxdt).*((-1).*utdt+uxdx)+(-1).*uydt.*uydx+uy.^2.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-2).*utdt+uxdx+uydy)))+ux.^3.*uy.*(utdy.*(utdt+(-1).*uxdx)+(-1).*utdx.*uxdy+2.*uxdt.*(uxdy+uydx)+(-1).*uydt.*(utdt+(-2).*uxdx+uydy)+uy.^2.*((-1).*(utdx+(-1).*uxdt).*(uxdy+uydx)+(-1).*(utdy+(-1).*uydt).*((-2).*utdt+uxdx+uydy)))+ux.^2.*(uxdt.*((-1).*utdt+uxdx)+(-1).*uydt.*uydx+uy.^4.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+uy.^2.*(utdx.*(utdt+(-1).*uxdx)+(-1).*uxdy.*(utdy+(-2).*uydt)+uxdt.*((-3).*utdt+uxdx+2.*uydy)))));
  t3dd(:,:,1,3) = -(1/2).*(ut.*(ux.^4.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2).*uy+2.*ux.^3.*(((-1).*utdx+uxdt).*uydt+((-1).*utdt+uxdx).*uydx+uy.^2.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+ux.^2.*uy.*(utdt.^2+(-1).*uxdx.^2+(-2).*utdy.*uydt+uydt.^2+2.*uxdy.*uydx+uydx.^2+2.*((-1).*utdt+uxdx).*uydy+uy.^2.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy)))+uy.*((-1).*uxdt.^2+uydt.^2+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+uy.^2.*((-1).*uxdt.^2+(-1).*uxdy.^2+2.*uydt.*((-1).*utdy+uydt)+((-1).*utdt+uydy).^2))+2.*ux.*(uxdt.*uydt+uy.^4.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+uy.^2.*((-1).*(utdx+(-2).*uxdt).*uydt+uydx.*((-1).*utdt+uydy)+uxdy.*((-1).*uxdx+uydy))))+2.*((-1).*ux.^5.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx).*uy+ux.^4.*(((-1).*utdt+uxdx).*uydt+((-1).*utdx+uxdt).*uydx+uy.^2.*((-1).*((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(-1).*(utdx+(-1).*uxdt).*(uxdy+uydx)))+(-1).*uy.^2.*(1+uy.^2).*(uxdt.*uxdy+(-1).*((-1).*utdy.*uy.^2+(1+uy.^2).*uydt).*((-1).*utdt+uydy))+ux.^3.*uy.*((-1).*utdy.*uydx+uydt.*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+uydy)+uy.^2.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-2).*utdt+uxdx+uydy)))+ux.^2.*(((-1).*utdt+uxdx).*uydt+uxdt.*uydx+uy.^4.*((-1).*(utdx+(-1).*uxdt).*(uxdy+uydx)+(-1).*(utdy+(-1).*uydt).*((-2).*utdt+uxdx+uydy))+uy.^2.*((-1).*(utdx+(-2).*uxdt).*uydx+utdy.*(utdt+(-1).*uydy)+uydt.*((-3).*utdt+2.*uxdx+uydy)))+ux.*uy.*(uydt.*(uxdy+uydx)+uxdt.*((-1).*uxdx+uydy)+uy.^4.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+uy.^2.*((-1).*utdy.*uydx+2.*uydt.*(uxdy+uydx)+utdx.*(utdt+(-1).*uydy)+uxdt.*((-1).*utdt+(-1).*uxdx+2.*uydy)))));
  t3dd(:,:,2,2) = (1/2).*(uxdt.^2+ux.^6.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2)+(-1).*uydt.^2+ux.^5.*((-2).*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+2.*uy.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+2.*ut.*uy.*(uxdt.*uxdy+uydt.*(utdt+(-1).*uydy))+2.*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+uy.^2.*(uxdt.^2+uxdy.^2+2.*(utdy+(-1).*uydt).*uydt+(-1).*((-1).*utdt+uydy).^2)+(-1).*uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+ux.^4.*((utdx+(-3).*uxdt).*(utdx+(-1).*uxdt)+2.*((-1).*utdt+uxdx).^2+(-1).*uydt.^2+(-1).*uydx.^2+(-2).*ut.*uy.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+uy.^2.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy)))+2.*ux.*(ut.*(uxdt.*((-1).*utdt+uxdx)+(-1).*uydt.*uydx)+uy.*((-1).*utdy.*uxdt+((-1).*utdt+uxdx).*uxdy+(utdx+uxdt).*uydt+uydx.*(utdt+(-1).*uydy))+ut.*uy.^2.*((-1).*(utdy+(-1).*uydt).*(uxdy+(-1).*uydx)+(utdx+uxdt).*((-1).*utdt+uydy))+uy.^3.*((-1).*(utdx+uxdt).*(utdy+(-1).*uydt)+(uxdy+(-1).*uydx).*((-1).*utdt+uydy)))+2.*ux.^3.*(ut.*((-1).*(utdx+(-2).*uxdt).*((-1).*utdt+uxdx)+(-1).*uydt.*uydx)+uy.*(utdy.*(utdx+(-2).*uxdt)+2.*((-1).*utdt+uxdx).*uxdy+2.*uxdt.*uydt+uydx.*(uxdx+(-1).*uydy))+(-1).*ut.*uy.^2.*((utdy+(-1).*uydt).*(uxdy+uydx)+(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+uy.^3.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy)))+ux.^2.*((-2).*utdx.*uxdt+3.*uxdt.^2+((-1).*utdt+uxdx).^2+(-2).*uydt.^2+(-1).*uydx.^2+2.*ut.*uy.*(utdy.*(utdt+(-1).*uxdx)+utdx.*((-1).*uxdy+uydx)+uxdt.*(2.*uxdy+uydx)+uydt.*(uxdx+(-1).*uydy))+(-2).*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+uy.^2.*((-1).*utdx.^2+utdy.^2+(-2).*utdx.*uxdt+2.*uxdt.^2+2.*uxdy.^2+(-1).*uydt.^2+2.*uxdy.*uydx+(-1).*uydx.^2+(-1).*((-1).*utdt+uydy).*(utdt+(-2).*uxdx+uydy))));
  t3dd(:,:,2,3) = (1/2).*(ux.^5.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2).*uy+2.*(1+uy.^2).*(ut.*uxdy.*uy+uxdt.*(1+uy.^2)).*uydt+2.*ux.^4.*((-1).*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx).*uy+((-1).*utdx+uxdt).*uydt+((-1).*utdt+uxdx).*uydx+uy.^2.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+(-2).*uxdt.*uy.*(1+uy.^2).*(utdy.*uy+ut.*(utdt+(-1).*uydy))+(-2).*uxdy.*uy.^2.*(ut.*utdy.*uy+(-1).*(1+uy.^2).*((-1).*utdt+uydy))+ux.*(2.*ut.*(((-1).*utdt+uxdx).*uydt+uxdt.*uydx)+(-2).*ut.*uy.^4.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+uy.*((-2).*utdx.*uxdt+uxdt.^2+(-2).*utdy.*uydt+uydt.^2+2.*uxdy.*uydx+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy))+uy.^5.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+2.*ut.*uy.^2.*((-1).*utdx.*uxdy+uxdt.*(uxdy+uydx)+(-1).*(utdy+(-1).*uydt).*((-2).*utdt+uxdx+uydy))+uy.^3.*((-2).*utdx.*uxdt+uxdt.^2+uxdy.^2+2.*(utdy+(-1).*uydt).^2+2.*uxdy.*uydx+((-1).*utdt+uydy).*((-3).*utdt+2.*uxdx+uydy)))+ux.^3.*(2.*ut.*(((-1).*utdt+uxdx).*uydt+((-1).*utdx+uxdt).*uydx)+(-2).*ut.*uy.^2.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+uy.^3.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy))+uy.*(2.*(utdx+(-1).*uxdt).^2+(-2).*utdy.*uydt+uydt.^2+2.*uxdy.*uydx+uydx.^2+((-1).*utdt+uxdx).*((-3).*utdt+uxdx+2.*uydy)))+2.*ux.^2.*((-1).*(utdx+(-2).*uxdt).*uydt+((-1).*utdt+uxdx).*uydx+uy.^4.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+uy.^2.*(2.*utdy.*(utdx+(-1).*uxdt)+((-2).*utdx+3.*uxdt).*uydt+(uxdy+uydx).*((-2).*utdt+uxdx+uydy))+ut.*((uy+uy.^3).*uydt.*(uxdy+uydx)+uy.*((-1).*utdy.*uydx+(-1).*(utdx+(-1).*uxdt).*((-2).*utdt+uxdx+uydy)+(-1).*uy.^2.*(utdy.*(uxdy+uydx)+(utdx+(-1).*uxdt).*((-1).*utdt+uydy))))));
  t3dd(:,:,3,3) = (1/2).*(ux.^4.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2).*((-1)+uy.^2)+ux.^3.*(2.*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+(-2).*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx).*uy.^2+(-2).*uy.*((utdx+(-1).*uxdt).*(utdy+uydt)+((-1).*utdt+uxdx).*(uxdy+(-1).*uydx))+2.*uy.^3.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+ux.^2.*(2.*(utdx+(-1).*uxdt).*uxdt+(-1).*((-1).*utdt+uxdx).^2+uydt.^2+2.*ut.*uy.*(((-1).*utdt+uxdx).*(utdy+uydt)+(utdx+(-1).*uxdt).*(uxdy+(-1).*uydx))+uydx.^2+(-2).*ut.*uy.^3.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+(-1).*uy.^2.*((-1).*utdx.^2+utdy.^2+uxdt.^2+uxdy.^2+2.*utdy.*uydt+(-2).*uydt.^2+(-2).*uxdy.*uydx+(-2).*uydx.^2+((-1).*utdt+uxdx).*(utdt+uxdx+(-2).*uydy))+uy.^4.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy)))+(1+uy.^2).*((-1).*uxdt.^2+uydt.^2+(-2).*ut.*uy.*(uxdt.*uxdy+uydt.*(utdt+(-1).*uydy))+(-2).*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+uy.^2.*((-1).*uxdt.^2+(-1).*uxdy.^2+2.*uydt.*((-1).*utdy+uydt)+((-1).*utdt+uydy).^2))+2.*ux.*(uy.*(1+uy.^2).*((utdt+(-1).*uxdx).*uxdy+(-1).*utdx.*uydt+uxdt.*(utdy+uydt)+uydx.*((-1).*utdt+uydy)+uy.^2.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy)))+ut.*(uxdt.*(utdt+(-1).*uxdx)+uydt.*uydx+uy.^4.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+uy.^2.*(utdy.*(uxdy+(-1).*uydx)+uydt.*(uxdy+2.*uydx)+utdx.*(utdt+(-1).*uydy)+uxdt.*((-1).*uxdx+uydy)))));

  t4dd(:,:,1,1) = (1/2).*(ux.^6.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2)+ux.^5.*((-2).*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+2.*uy.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+uy.^2.*((-1).*uxdx.^2+(-1).*uxdy.^2+uydx.^2+uydy.^2+(-2).*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+(-2).*ut.*uy.*(uxdt.*uxdy+utdx.*uydx+(utdy+(-1).*uydt).*uydy)+uy.^2.*(utdx.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(utdy+(-1).*uydt).^2+uydx.^2+2.*uydy.*((-1).*utdt+uydy))+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2))+ux.^3.*((-2).*ut.*((utdx+(-1).*uxdt).*uxdx+utdy.*uxdy+uydt.*uydx)+(-2).*ut.*uy.^2.*((utdy+(-1).*uydt).*(uxdy+uydx)+(utdx+(-1).*uxdt).*((-2).*utdt+uxdx+uydy))+2.*uy.^3.*(2.*(utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-2).*utdt+uxdx+uydy))+2.*uy.*((-1).*utdy.*uxdt+(-1).*(utdx+(-2).*uxdt).*uydt+(-1).*uydx.*(utdt+(-3).*uxdx+uydy)+uxdy.*((-1).*utdt+uxdx+uydy)))+ux.^2.*(uxdx.^2+uxdy.^2+(-1).*uydx.^2+(-1).*uydy.^2+(-2).*ut.*uy.*(utdx.*uydx+(-1).*uxdt.*(uxdy+2.*uydx)+utdy.*uydy+uydt.*((-2).*uxdx+uydy))+2.*uy.^2.*(utdx.*(utdx+(-1).*uxdt)+utdy.*(utdy+(-1).*uydt)+(uxdy+uydx).^2+2.*uxdx.*uydy+(-1).*utdt.*(uxdx+uydy))+2.*ut.*uy.^3.*((-1).*(utdx+(-1).*uxdt).*(uxdy+uydx)+(-1).*(utdy+(-1).*uydt).*((-2).*utdt+uxdx+uydy))+uy.^4.*((utdx+(-1).*uxdt).^2+2.*(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+((-1).*utdt+uydy).*((-3).*utdt+2.*uxdx+uydy)))+ux.^4.*(utdx.^2+utdy.^2+(-2).*utdx.*uxdt+uxdt.^2+(-2).*utdt.*uxdx+2.*uxdx.^2+uxdy.^2+(-1).*uydt.^2+(-1).*uydx.^2+(-2).*ut.*uy.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+uy.^2.*(2.*(utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+((-1).*utdt+uxdx).*((-3).*utdt+uxdx+2.*uydy)))+ux.*uy.*(4.*uxdx.*uydx+4.*uxdy.*uydy+(-2).*ut.*uy.*((utdx+uxdt).*uxdx+utdy.*uxdy+(-1).*uydt.*(2.*uxdy+uydx)+(-2).*uxdt.*uydy)+(-2).*ut.*uy.^3.*((utdy+(-1).*uydt).*(uxdy+uydx)+(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+2.*uy.^4.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+2.*uy.^2.*((-1).*utdy.*uxdt+(-1).*(utdx+(-2).*uxdt).*uydt+uydx.*((-1).*utdt+uxdx+uydy)+uxdy.*((-1).*utdt+(-1).*uxdx+3.*uydy))));
  t4dd(:,:,1,2) = -(1/2).*(ut.*(ux.^5.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2)+ux.^3.*(utdy.^2+(utdx+(-1).*uxdt).^2+2.*uxdx.*((-1).*utdt+uxdx)+uxdy.^2+(-1).*uydt.^2+(-1).*uydx.^2)+2.*ux.^4.*uy.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx))+2.*uy.*(uxdx.*uydx+uxdy.*uydy)+ux.*(uxdx.^2+uxdy.^2+(-1).*uydx.^2+(-1).*uydy.^2)+ux.*uy.^2.*(utdy.^2+(utdx+(-1).*uxdt).^2+(-1).*uydt.^2+(uxdy+uydx).^2+2.*uxdx.*((-1).*utdt+uydy))+ux.^3.*uy.^2.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy))+2.*uy.^3.*(uxdt.*((-1).*utdy+uydt)+uxdy.*((-1).*utdt+uydy))+2.*ux.^2.*uy.^3.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+ux.*uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+2.*ux.^2.*uy.*(uxdt.*((-1).*utdy+uydt)+(-1).*uydx.*((-2).*uxdx+uydy)+uxdy.*((-1).*utdt+uxdx+uydy)))+2.*((-1).*ux.^6.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+ux.^5.*uy.*((-1).*((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(-1).*(utdx+(-1).*uxdt).*(uxdy+uydx))+uy.^2.*((-1).*utdx.*uxdx+(-1).*utdy.*uxdy.*(1+uy.^2)+uxdy.*(1+uy.^2).*uydt+uxdt.*((-1).*utdt.*uy.^2+(1+uy.^2).*uydy))+ux.^3.*uy.*((-1).*(utdx.*(1+uy.^2)+(-1).*uxdt.*(2+uy.^2)).*(uxdy+uydx)+uydt.*((-2).*utdt.*uy.^2+uxdx.*(2+uy.^2)+((-1)+uy.^2).*uydy)+utdy.*(utdt.*(1+2.*uy.^2)+(-1).*(1+uy.^2).*(uxdx+uydy)))+ux.^4.*(uydt.*(uxdy.*uy.^2+((-1)+uy.^2).*uydx)+(-1).*utdy.*(uxdy+uy.^2.*(uxdy+uydx))+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+2.*uxdx+uy.^2.*((-2).*utdt+uxdx+uydy)))+ux.*uy.*(uxdt.*(uxdy+uydx)+uydt.*(uxdx+(-1).*uydy)+(-1).*uy.^4.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+uy.^2.*(uxdx.*uydt+(-1).*(utdx+(-1).*uxdt).*(uxdy+uydx)+(-1).*utdy.*((-1).*utdt+uxdx+uydy)))+ux.^2.*((-1).*utdx.*uxdx+uxdt.*uxdx+(-1).*utdy.*uxdy+(-1).*uydt.*uydx+uy.^4.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+uy.^2.*(utdx.*(utdt+(-2).*uxdx)+2.*uxdy.*((-1).*utdy+uydt)+uxdt.*((-2).*utdt+uxdx+2.*uydy)))));
  t4dd(:,:,1,3) = -(1/2).*(ut.*(ux.^4.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2).*uy+2.*ux.^3.*(((-1).*utdx+uxdt).*uydt+((-1).*utdt+uxdx).*uydx)+2.*ux.^3.*uy.^2.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx))+ux.^2.*uy.*(utdx.^2+(-1).*uxdt.^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*uydy)+2.*ux.*(uxdx.*uydx+uxdy.*uydy)+uy.*((-1).*uxdx.^2+(-1).*uxdy.^2+uydx.^2+uydy.^2)+ux.^2.*uy.^3.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy))+2.*ux.*uy.^4.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+uy.^3.*(utdx.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(utdy+(-1).*uydt).^2+uydx.^2+2.*uydy.*((-1).*utdt+uydy))+uy.^5.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+2.*ux.*uy.^2.*(((-1).*utdx+uxdt).*uydt+uydx.*((-1).*utdt+uxdx+uydy)+uxdy.*((-1).*uxdx+2.*uydy)))+2.*((-1).*ux.^5.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx).*uy+ux.^4.*(((-1).*utdt+uxdx).*uydt+((-1).*utdx+uxdt).*uydx+uy.^2.*((-1).*((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(-1).*(utdx+(-1).*uxdt).*(uxdy+uydx)))+(-1).*uy.^2.*(1+uy.^2).*(uxdt.*uxdy+utdx.*uydx+(utdy+(-1).*uydt).*((-1).*utdt.*uy.^2+(1+uy.^2).*uydy))+ux.^3.*uy.*(utdx.*(utdt+(-1).*uxdx)+(-1).*utdy.*(1+uy.^2).*(uxdy+uydx)+uydt.*(uxdy+uydx)+(-1).*utdx.*uydy+uxdt.*uydy+uy.^2.*(uydt.*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-2).*utdt+uxdx+uydy)))+ux.^2.*((-1).*(utdx+(-1).*uxdt).*(uxdy.*uy.^4+(1+uy.^2).^2.*uydx)+(-1).*utdy.*(uxdx.*uy.^4+(-1).*utdt.*uy.^2.*(1+2.*uy.^2)+(1+uy.^2).^2.*uydy)+(1+uy.^2).*uydt.*(uxdx+uy.^2.*((-2).*utdt+uxdx+uydy)))+ux.*uy.*(uydt.*(uxdy+uydx)+uxdt.*((-1).*uxdx+uydy)+uy.^4.*((-1).*(utdy+(-1).*uydt).*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+uy.^2.*((-1).*(utdy+(-2).*uydt).*(uxdy+uydx)+(-1).*utdx.*((-1).*utdt+uxdx+uydy)+uxdt.*((-1).*uxdx+2.*uydy)))));
  t4dd(:,:,2,2) = (1/2).*(uxdx.^2+ux.^6.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2)+uxdy.^2+(-1).*uydx.^2+ux.^5.*((-2).*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+2.*uy.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+(-1).*uydy.^2+2.*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+2.*ut.*uy.*(uxdt.*uxdy+utdx.*uydx+(utdy+(-1).*uydt).*uydy)+(-1).*uy.^2.*(utdx.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(utdy+(-1).*uydt).^2+uydx.^2+2.*uydy.*((-1).*utdt+uydy))+(-1).*uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+ux.^4.*(utdt.^2+2.*utdx.^2+utdy.^2+(-4).*utdx.*uxdt+2.*uxdt.^2+(-4).*utdt.*uxdx+3.*uxdx.^2+uxdy.^2+(-1).*uydt.^2+(-1).*uydx.^2+(-2).*ut.*uy.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+uy.^2.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy)))+ux.^2.*(utdx.^2+utdy.^2+(-2).*utdx.*uxdt+uxdt.^2+(-2).*utdt.*uxdx+3.*uxdx.^2+2.*uxdy.^2+(-1).*uydt.^2+(-2).*uydx.^2+(-1).*uydy.^2+(-2).*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+2.*uy.^2.*(utdy.^2+(-1).*utdx.*uxdt+uxdt.^2+(-1).*uydt.^2+uxdy.*(uxdy+uydx)+uxdx.*((-1).*utdt+uydy))+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+2.*ut.*uy.*((-1).*utdx.*uxdy+uxdt.*(2.*uxdy+uydx)+uydt.*(utdt+uxdx+(-1).*uydy)+(-1).*utdy.*((-1).*utdt+uxdx+uydy)))+ux.*((-2).*ut.*((utdx+(-1).*uxdt).*uxdx+utdy.*uxdy+uydt.*uydx)+(-2).*ut.*uy.^2.*((utdy+(-1).*uydt).*(uxdy+(-1).*uydx)+(-1).*(utdx+uxdt).*((-1).*utdt+uydy))+2.*uy.^3.*((-1).*(utdx+uxdt).*(utdy+(-1).*uydt)+(uxdy+(-1).*uydx).*((-1).*utdt+uydy))+2.*uy.*((-1).*utdy.*uxdt+utdx.*uydt+uydx.*(utdt+uxdx+(-1).*uydy)+uxdy.*((-1).*utdt+uxdx+uydy)))+ux.^3.*((-2).*ut.*((utdx+(-1).*uxdt).*((-1).*utdt+2.*uxdx)+utdy.*uxdy+uydt.*uydx)+(-2).*ut.*uy.^2.*((utdy+(-1).*uydt).*(uxdy+uydx)+(utdx+(-1).*uxdt).*((-1).*utdt+uydy))+2.*uy.^3.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+2.*uy.*(utdy.*(utdx+(-2).*uxdt)+uxdt.*uydt+(-1).*uydx.*((-2).*uxdx+uydy)+uxdy.*((-2).*utdt+2.*uxdx+uydy))));
  t4dd(:,:,2,3) = (1/2).*(ux.^5.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2).*uy+2.*ux.^4.*((-1).*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx).*uy+((-1).*utdx+uxdt).*uydt+((-1).*utdt+uxdx).*uydx+uy.^2.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+ux.^3.*(2.*ut.*(((-1).*utdt+uxdx).*uydt+((-1).*utdx+uxdt).*uydx)+(-2).*ut.*uy.^2.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+uy.^3.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy))+uy.*(3.*utdx.^2+(-4).*utdx.*uxdt+uxdt.^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uxdx+uydy)))+ux.*((-2).*ut.*uy.^4.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+2.*ut.*(uxdx.*uydt+((-1).*utdx+uxdt).*uydx+(-1).*utdy.*uydy)+uy.^5.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2)+uy.*(2.*utdx.*(utdx+(-1).*uxdt)+2.*utdy.*(utdy+(-1).*uydt)+(uxdy+uydx).^2+(uxdx+uydy).*((-2).*utdt+uxdx+uydy))+uy.^3.*(3.*utdy.^2+(utdx+(-1).*uxdt).^2+(-4).*utdy.*uydt+uydt.^2+(uxdy+uydx).^2+2.*((-1).*utdt+uydy).*((-1).*utdt+uxdx+uydy))+2.*ut.*uy.^2.*((-1).*(utdx+(-1).*uxdt).*(uxdy+uydx)+uydt.*((-1).*utdt+uxdx+uydy)+(-1).*utdy.*((-2).*utdt+uxdx+2.*uydy)))+2.*(ut.*uy.*((-1).*utdx.*uxdx+(-1).*utdy.*uxdy.*(1+uy.^2)+uxdy.*(1+uy.^2).*uydt+uxdt.*((-1).*utdt.*uy.^2+(1+uy.^2).*uydy))+(1+uy.^2).*(uxdt.*uy.^2.*((-1).*utdy+uydt)+uxdx.*uydx+uxdy.*((-1).*utdt.*uy.^2+(1+uy.^2).*uydy)))+2.*ux.^2.*((-1).*utdx.*uydt+uxdt.*uydt+(-1).*utdt.*uydx+2.*uxdx.*uydx+uxdy.*uydy+uy.^4.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+(uxdy+uydx).*((-1).*utdt+uydy))+uy.^2.*(2.*(utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+uydx.*((-2).*utdt+2.*uxdx+uydy)+uxdy.*((-2).*utdt+uxdx+2.*uydy))+ut.*((uy+uy.^3).*uydt.*(uxdy+uydx)+uy.*((-1).*utdy.*(1+uy.^2).*(uxdy+uydx)+uxdt.*(uxdx+(-1).*utdt.*(1+uy.^2)+(1+uy.^2).*uydy)+(-1).*utdx.*(2.*uxdx+(-1).*utdt.*(2+uy.^2)+(1+uy.^2).*uydy)))));
  t4dd(:,:,3,3) = (1/2).*(ux.^4.*((utdx+(-1).*uxdt).^2+((-1).*utdt+uxdx).^2).*((-1)+uy.^2)+ux.^3.*(2.*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx)+(-2).*ut.*(utdx+(-1).*uxdt).*((-1).*utdt+uxdx).*uy.^2+(-2).*uy.*((utdx+(-1).*uxdt).*(utdy+uydt)+((-1).*utdt+uxdx).*(uxdy+(-1).*uydx))+2.*uy.^3.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+((-1).*utdt+uxdx).*(uxdy+uydx)))+ux.^2.*((-1).*utdx.^2+(-1).*utdy.^2+2.*utdx.*uxdt+(-1).*uxdt.^2+2.*utdt.*uxdx+(-2).*uxdx.^2+(-1).*uxdy.^2+uydt.^2+2.*ut.*uy.*(((-1).*utdt+uxdx).*(utdy+uydt)+(utdx+(-1).*uxdt).*(uxdy+(-1).*uydx))+uydx.^2+(-2).*ut.*uy.^3.*(((-1).*utdt+uxdx).*(utdy+(-1).*uydt)+(utdx+(-1).*uxdt).*(uxdy+uydx))+2.*uy.^2.*(utdx.^2+(-1).*uxdt.^2+(-1).*utdy.*uydt+uydt.^2+uydx.*(uxdy+uydx)+((-1).*utdt+uxdx).*uydy)+uy.^4.*((utdx+(-1).*uxdt).^2+(utdy+(-1).*uydt).^2+(uxdy+uydx).^2+2.*((-1).*utdt+uxdx).*((-1).*utdt+uydy)))+(1+uy.^2).*((-1).*uxdx.^2+(-1).*uxdy.^2+uydx.^2+uydy.^2+(-2).*ut.*uy.^3.*(utdy+(-1).*uydt).*((-1).*utdt+uydy)+(-2).*ut.*uy.*(uxdt.*uxdy+utdx.*uydx+(utdy+(-1).*uydt).*uydy)+uy.^2.*(utdx.^2+(-1).*uxdt.^2+(-1).*uxdy.^2+(utdy+(-1).*uydt).^2+uydx.^2+2.*uydy.*((-1).*utdt+uydy))+uy.^4.*((utdy+(-1).*uydt).^2+((-1).*utdt+uydy).^2))+2.*ux.*(ut.*((utdx+(-1).*uxdt).*uxdx+uydt.*uydx+(-1).*utdy.*(1+uy.^2).*(uxdy.*((-1)+uy.^2)+uy.^2.*uydx)+uy.^2.*((-1).*(utdx+uxdt).*((-1).*utdt+uxdx)+uydt.*(uxdy+2.*uydx)+((-1).*utdx+uxdt).*uydy)+uy.^4.*(uydt.*(uxdy+uydx)+(-1).*(utdx+(-1).*uxdt).*((-1).*utdt+uydy)))+uy.*(1+uy.^2).*(utdy.*uxdt+(-1).*utdx.*uydt+uydx.*((-1).*utdt+uxdx+uydy)+uxdy.*((-1).*uxdx+(-1).*utdt.*((-1)+uy.^2)+(1+uy.^2).*uydy)+uy.^2.*((utdx+(-1).*uxdt).*(utdy+(-1).*uydt)+uydx.*((-1).*utdt+uydy)))));

  t5dd(:,:,1,1) = (1/2).*(ux.^4.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+ux.^3.*(2.*ut.*((-1).*utdt.*utdx+uxdt.*uxdx+uydt.*uydx)+2.*uy.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy))+uy.^2.*(utdx.^2+(-1).*utdy.^2+(-1).*uxdx.^2+uxdy.^2+(-1).*uydx.^2+uydy.^2+2.*ut.*uy.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+uy.^2.*((-1).*utdt.^2+(-1).*utdy.^2+uxdt.^2+uxdy.^2+uydt.^2+uydy.^2))+ux.^2.*((-1).*utdx.^2+utdy.^2+uxdx.^2+(-1).*uxdy.^2+uydx.^2+(-1).*uydy.^2+2.*ut.*uy.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+uy.^2.*((-1).*utdx.^2+(-1).*utdy.^2+uxdx.^2+uxdy.^2+2.*((-1).*utdt.^2+uxdt.^2+uydt.^2)+uydx.^2+uydy.^2))+2.*ux.*uy.*((-1).*utdx.*utdy.*(2+uy.^2)+uxdx.*uxdy.*(2+uy.^2)+2.*uydx.*uydy+uy.*(ut.*((-1).*utdt.*utdx+uxdt.*uxdx)+uydx.*(ut.*uydt+uy.*uydy))));
  t5dd(:,:,1,2) = -(1/2).*((-2).*(ux.^2+uy.^2).*(utdt.*utdx+(-1).*uxdt.*uxdx+(-1).*uydt.*uydx+ux.^2.*(utdt.*utdx+(-1).*uxdt.*uxdx+(-1).*uydt.*uydx)+ux.*uy.*(utdt.*utdy+(-1).*uxdt.*uxdy+(-1).*uydt.*uydy))+ut.*(ux.^3.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+2.*uy.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy)+2.*ux.^2.*uy.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy)+ux.*((-1).*utdx.^2+uxdx.^2+(-1).*utdy.^2.*((-1)+uy.^2)+uxdy.^2.*((-1)+uy.^2)+uydx.^2+(-1).*uydy.^2+uy.^2.*((-1).*utdt.^2+uxdt.^2+uydt.^2+uydy.^2))));
  t5dd(:,:,1,3) = -(1/2).*((-2).*(ux.^2+uy.^2).*(utdt.*utdy.*(1+uy.^2)+ux.*uy.*(utdt.*utdx+(-1).*uxdt.*uxdx+(-1).*uydt.*uydx)+(-1).*(1+uy.^2).*(uxdt.*uxdy+uydt.*uydy))+ut.*(2.*ux.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy)+2.*ux.*uy.^2.*((-1).*utdx.*utdy+uxdx.*uxdy+uydx.*uydy)+uy.^3.*((-1).*utdt.^2+(-1).*utdy.^2+uxdt.^2+uxdy.^2+uydt.^2+uydy.^2)+uy.*((-1).*utdy.^2+(-1).*utdx.^2.*((-1)+ux.^2)+((-1)+ux.^2).*uxdx.^2+uxdy.^2+(-1).*uydx.^2+ux.^2.*((-1).*utdt.^2+uxdt.^2+uydt.^2+uydx.^2)+uydy.^2)));
  t5dd(:,:,2,2) = (1/2).*((-1).*utdx.^2+uxdx.^2+(-1).*uxdy.^2+2.*utdy.*((-1).*ut.*utdt.*((-1)+ux.^2)+(-1).*utdx.*(ux+ux.^3)).*uy+utdy.^2.*(1+uy.^2+(-1).*ux.^2.*((-1)+uy.^2))+uydx.^2+ux.^4.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+(-1).*uydy.^2+(-2).*ut.*uy.*(uxdt.*uxdy+uydt.*uydy)+(-1).*uy.^2.*((-1).*utdt.^2+uxdt.^2+uxdy.^2+uydt.^2+uydy.^2)+ux.^2.*((-2).*utdx.^2+2.*uxdx.^2+2.*ut.*uxdt.*uxdy.*uy+uxdy.^2.*((-1)+uy.^2)+(-1).*(1+uy.^2).*(utdt.^2+(-1).*uxdt.^2+(-1).*uydt.^2)+2.*uydx.^2+2.*ut.*uy.*uydt.*uydy+((-1)+uy.^2).*uydy.^2)+ux.*(2.*ut.*((-1).*utdt.*utdx+uxdt.*uxdx+uydt.*uydx)+2.*uy.*(uxdx.*uxdy+uydx.*uydy))+ux.^3.*(2.*ut.*((-1).*utdt.*utdx+uxdt.*uxdx+uydt.*uydx)+2.*uy.*(uxdx.*uxdy+uydx.*uydy)));
  t5dd(:,:,2,3) = (1/2).*(ux.^3.*uy.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+ux.*(2.*ut.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+2.*ut.*uy.^2.*((-1).*utdt.*utdy+uxdt.*uxdy+uydt.*uydy)+uy.^3.*((-1).*utdt.^2+(-1).*utdy.^2+uxdt.^2+uxdy.^2+uydt.^2+uydy.^2)+uy.*((-1).*utdx.^2+(-1).*utdy.^2+uxdx.^2+uxdy.^2+2.*((-1).*utdt.^2+uxdt.^2+uydt.^2)+uydx.^2+uydy.^2))+2.*((-1).*utdx.*utdy.*(1+uy.^2)+uxdx.*uxdy.*(1+uy.^2)+uydx.*uydy+uy.*(ut.*((-1).*utdt.*utdx+uxdt.*uxdx)+uydx.*(ut.*uydt+uy.*uydy)))+2.*ux.^2.*((-1).*utdx.*utdy.*(1+uy.^2)+uxdx.*uxdy.*(1+uy.^2)+uydx.*uydy+uy.*(ut.*((-1).*utdt.*utdx+uxdt.*uxdx)+uydx.*(ut.*uydt+uy.*uydy))));
  t5dd(:,:,3,3) = (1/2).*(utdx.^2.*(1+ux.^2)+(-1).*uxdx.^2+(-2).*utdy.*(ut.*utdt+utdx.*ux).*uy.*(1+uy.^2)+2.*(ut.*uxdt+ux.*uxdx).*uxdy.*uy.*(1+uy.^2)+(-1).*utdy.^2.*(1+uy.^2).^2+uxdy.^2.*(1+uy.^2).^2+(-1).*uydx.^2+ut.*ux.*(2.*utdt.*utdx+(-2).*uxdt.*uxdx+(-2).*uydt.*uydx)+2.*ut.*ux.*uy.^2.*((-1).*utdt.*utdx+uxdt.*uxdx+uydt.*uydx)+(-1).*ux.^2.*((-1).*utdt.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+ux.^2.*uy.^2.*((-1).*utdt.^2+(-1).*utdx.^2+uxdt.^2+uxdx.^2+uydt.^2+uydx.^2)+2.*ut.*uy.*uydt.*uydy+2.*ut.*uy.^3.*uydt.*uydy+2.*ux.*uy.*uydx.*uydy+2.*ux.*uy.^3.*uydx.*uydy+uydy.^2+uy.^4.*((-1).*utdt.^2+uxdt.^2+uydt.^2+uydy.^2)+uy.^2.*((-1).*utdt.^2+utdx.^2+uxdt.^2+(-1).*uxdx.^2+uydt.^2+(-1).*uydx.^2+2.*uydy.^2));


  
  % ****************************************************
  % ****************************************************

  ud = zeros(Nx, Ny, 3);
  ud(:,:,1) = -ut;
  ud(:,:,2) =  ux;
  ud(:,:,3) =  uy;

  uudd = zeros(Nx, Ny, 3, 3);
  uudd(:,:,1,1) =  ut.*ut;
  uudd(:,:,1,2) = -ut.*ux;
  uudd(:,:,1,3) = -ut.*uy;
  uudd(:,:,2,2) =  ux.*ux;
  uudd(:,:,2,3) =  ux.*uy;
  uudd(:,:,3,3) =  uy.*uy;

  Pdd = uudd;
  Pdd(:,:,1,1) = Pdd(:,:,1,1) - 1;
  Pdd(:,:,2,2) = Pdd(:,:,2,2) + 1;
  Pdd(:,:,3,3) = Pdd(:,:,3,3) + 1;






  % ****************************************************
  % ****************************************************
  %  flip signs
  ut = -ut;

  utut = ut.*ut;
  utux = ut.*ux;
  utuy = ut.*uy;
  uxux = ux.*ux;
  uxuy = ux.*uy;
  uyuy = uy.*uy;

  % ****************************************************
  % ****************************************************

  g = zeros(params.Nx, params.Ny, Nz, 4, 4);
  if true
  for iz = 1:Nz
    z = params.Cz(iz);
    % *************************************************************
    % zeroth order: 
    % *************************************************************
    % compute g/z^2 at zeroth order
    % compute drg/z^3

    % term numbers refer to eq (21) in Marks paper

    % term 1: -2 u_mu   dx^mu dr
    %  dr = -1/z^2 dz
    % term 1:  2 u_mu / z^2  dx^mu dz
    g(:,:,iz,1:3,4) = reshape(ud, [Nx,Ny,1,3,1]);
    g(:,:,iz,  4,4) = 0;

    % term 2 only
    % -1/z^2  (1 - (z T)^3)

    f = 1-(z.*T).^3;
    % size(f)    = [Nx Ny]
    % size(uudd) = [Nx Ny 3 3]
    g(:,:,iz,1:3,1:3) = reshape(  repmat(-f, [1,1,3,3]) .* uudd  , [Nx,Ny,1,3,3]);

    % term 3 only
    % 1/z^2
    g(:,:,iz,1:3,1:3) = g(:,:,iz,1:3,1:3) + reshape( Pdd, [Nx,Ny,1,3,3]);
  end
  end



  % *************************************************************
  % first order: 
  % *************************************************************

  if params.metricOrder >= 1
    disp('first');
    for iz = 1:Nz
      z = params.Cz(iz);
      % remember g has 1/z^2 factored out!

      % ********************************************************
      % h terms 

      % nothing!

      % ********************************************************
      % k terms 

      if (true) % checks
        Fk1 = z.*delu; % with 1/z^4 extracted

        % k/r^2 u u  dx dx
        % k z^2 u u dx dx
        % 1/z^2 k' u u 
        g(:,:,iz,1:3,1:3) = g(:,:,iz,1:3,1:3) + reshape(  repmat(Fk1, [1,1,3,3]) .* uudd  , [Nx,Ny,1,3,3]);
      end

      % ********************************************************
      % j terms 

      if (true)
        % j_mu = 1/z^2 u^lam d_lam d_mu 
        % -z(j_mu u_nu + j_nu u_mu)
        % -1/z j' u
        % 1/z^2(-z j' u)

        Fj1 = -z;

        for i=1:3
          for j=i:3
            g(:,:,iz,i,j) = g(:,:,iz,i,j) + Fj1.*(j1d(:,:,i).*ud(:,:,j) + j1d(:,:,j).*ud(:,:,i));
          end
        end
      end

      % ********************************************************
      % a terms 

      if (true)
        %  br = r./T;
        F = (1./6).*(3.^(1./2).*pi+(-2).*3.^(1./2).*atan(3.^(-1./2).*(1+2.*T.^(-1).*z.^(-1)))+3.*log(1+T.*z+T.^2.*z.^2));

        Fa1 = 2.*F./T;

        g(:,:,iz,1:3,1:3) = g(:,:,iz,1:3,1:3) + reshape(  repmat(Fa1, [1,1,3,3]) .* a1dd,  [Nx,Ny,1,3,3]);
      end
    end
  end % end of the 0th and 1st order loop

  for i=1:4
    for j=i+1:4
      g(:,:,:,j,i) = g(:,:,:,i,j);
    end
  end

  detg1 = g(:,:,:,1,3).^2.*g(:,:,:,2,4).^2+(-2).*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,2,3)+g(:,:,:,1,4).^2.*g(:,:,:,2,3).^2+(-2).*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,3,4)+2.*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,2).*g(:,:,:,3,4)+(-2).*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,3).*g(:,:,:,3,4)+2.*g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,2).^2.*g(:,:,:,3,4).^2+(-1).*g(:,:,:,1,1).*g(:,:,:,2,2).*g(:,:,:,3,4).^2+2.*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,4).*g(:,:,:,3,3)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,4).^2.*g(:,:,:,3,3)+(-1).*g(:,:,:,1,4).^2.*g(:,:,:,2,2).*g(:,:,:,3,3);



  % *************************************************************
  % second order 
  % *************************************************************

  if (params.metricOrder >= 2 )
    %*********************************************
    % compute z-dependent functions first

    % initialize all tensor structures
    Fh2s6 = zeros(Nx, Ny, Nz);
    Fh2s7 = zeros(Nx, Ny, Nz);

    Fk2s3 = zeros(Nx, Ny, Nz);
    Fk2s4 = zeros(Nx, Ny, Nz);
    Fk2s5 = zeros(Nx, Ny, Nz);
    Fk2s6 = zeros(Nx, Ny, Nz);
    Fk2s7 = zeros(Nx, Ny, Nz);

    Fj2v3 = zeros(Nx, Ny, Nz);
    Fj2v4 = zeros(Nx, Ny, Nz);
    Fj2v5 = zeros(Nx, Ny, Nz);
    Fj2v6 = zeros(Nx, Ny, Nz);

    Fa2t2 = zeros(Nx, Ny, Nz); 
    Fa2t3 = zeros(Nx, Ny, Nz);
    Fa2t4 = zeros(Nx, Ny, Nz);
    Fa2t5 = zeros(Nx, Ny, Nz);

    for ix=1:Nx
      for iy=1:Ny
        z = params.Cz';

        Fh2s6(ix,iy,:) = h2s6(z, T(ix,iy));
        Fh2s7(ix,iy,:) = h2s7(z, T(ix,iy)); 

        Fk2s3(ix,iy,:) = k2s3(z, T(ix,iy));
        Fk2s4(ix,iy,:) = k2s4(z, T(ix,iy));
        Fk2s5(ix,iy,:) = k2s5(z, T(ix,iy));
        Fk2s6(ix,iy,:) = k2s6(z, T(ix,iy));
        Fk2s7(ix,iy,:) = k2s7(z, T(ix,iy));

        Fj2v3(ix,iy,:) = j2v3(z, T(ix,iy)); %a
        Fj2v4(ix,iy,:) = j2v4(z, T(ix,iy)); %b
        Fj2v5(ix,iy,:) = j2v5(z, T(ix,iy)); %c
        Fj2v6(ix,iy,:) = j2v6(z, T(ix,iy));

        Fa2t2(ix,iy,:) = a2t2(z, T(ix,iy)); %a
        Fa2t3(ix,iy,:) = Fa2t2(ix,iy,:);
        Fa2t4(ix,iy,:) = a2t4(z, T(ix,iy)); %b 
        Fa2t5(ix,iy,:) = a2t5(z, T(ix,iy)); %c
      end
    end

    % compute the actual metric
    for iz=1:Nz
      z = params.Cz(iz);

      % ********************************************************
      % h terms 

      % z^2 extracted (but has the factors of T inside)

      Fh2 = Fh2s6(:,:,iz) .* s6 + Fh2s7(:,:,iz) .* s7;
      % this has z^2 extracted

      % -2 b^2 h u_mu  dx^\nu dr
      %  dr = -1/z^2 dz!!!!
      % 1/z^2 2 b^2 h u_mu  dx^\nu dz
      % 1/z^2( 2 b^2 z^2 h' u_mu)  dx^\nu dz
      F = Fh2.*(z.^2)./(T.^2); % ok
      %                                              ok
      g(:,:,iz,1:3,4) = g(:,:,iz,1:3,4) + reshape(  repmat(F, [1,1,3]).*ud,    [Nx,Ny,1,3]);

      % -r^2 b^2 h P_mu_nu x^\mu x^\nu
      % -1/z^2 b^2 z^2 h' P_mu_nu x^\mu x^\nu
      % -b^2 h' P_mu_nu x^\mu x^\nu
      % 1/z^2(  -z ^2 b^2 h' P_mu_nu x^\mu x^\nu
      %  2 powers of z for overall h extraction
      % should be a minus
      F = -Fh2./(T.^2).*(z.^2);

      %                                                   ok
      g(:,:,iz,1:3,1:3) = g(:,:,iz,1:3,1:3) + reshape(  repmat(F, [1,1,3,3]) .* Pdd,  [Nx,Ny,1,3,3]);



      % ********************************************************
      % k terms 
      Fk2 = Fk2s3(:,:,iz).*s3 + Fk2s4(:,:,iz).*s4 + Fk2s5(:,:,iz).*s5 + Fk2s6(:,:,iz).*s6 + Fk2s7(:,:,iz).*s7; % 1/z^2 extracted

      % K + r/2 (1-4 b^3 r^3)h
      % K + 1/2/z (1-4/T^3/z^3) h
      % K + z/2 (1 - 4/T^3/z^3) h'
      % 1/z^2(K' + z^3/2(1-4/T^3/z^3) h'
      % 1/z^2(K' + ( z^3/2 - 2/T^3 ) h'
      Fk2 = Fk2 + ( (z.^3)./2 - 2./T.^3 ).*Fh2; % this has 1/z^2 extracted

      % k / b^2/r^2 u_mu u_nu term
      % 1/z^2 k'/ b^2  z^2
      % 1/z^2 ( z^2 k'/b^2 )
      F = Fk2.*(T.^2).*(z.^2);

      g(:,:,iz,1:3,1:3) = g(:,:,iz,1:3,1:3) + reshape(  repmat(F, [1,1,3,3]) .* uudd,  [Nx,Ny,1,3,3]);

      % ********************************************************
      % j terms 

      %  j2c with 1/z extracted!!

      j2d = zeros(Nx, Ny, 3);
      j2d(:,:,:) = ... 
        repmat(Fj2v3(:,:,iz), [1,1,3]) .* v3d + ...
        repmat(Fj2v4(:,:,iz), [1,1,3]) .* v4d + ...
        repmat(Fj2v5(:,:,iz), [1,1,3]) .* v5d + ...
        repmat(Fj2v6(:,:,iz), [1,1,3]) .* v6d;

      % -2/b/r j2 u_mu
      %  1/z from Fj, z from 1/r, have to pull out 1/z^2, so multiply by z^2
      % -2/b*z j2 
      % -2/b*z 1/z j2'
      % 1/z^2(-z^2 2/b j')
      F = -1.*z.^2.*T;

      for i=1:3
        for j=i:3
          g(:,:,iz,i,j) = g(:,:,iz,i,j) + F.*(j2d(:,:,i).*ud(:,:,j) + j2d(:,:,j).*ud(:,:,i));
        end
      end



      % ********************************************************
      % a terms 

      % checked block for typos 11/27
      a2dd = ...
        repmat(Fa2t2(:,:,iz), [1,1,3,3]).*t2dd + ...
        repmat(Fa2t3(:,:,iz), [1,1,3,3]).*t3dd + ...
        repmat(Fa2t4(:,:,iz), [1,1,3,3]).*t4dd + ...
        repmat(Fa2t5(:,:,iz), [1,1,3,3]).*t5dd;
        
        
      % b^2 r^2 a2
      F = 1./(T.^2);

      % checked block for typos
      g(:,:,iz,1:3,1:3) = g(:,:,iz,1:3,1:3) + reshape(  repmat(F, [1,1,3,3]) .* a2dd,  [Nx,Ny,1,3,3]);
    end
  end % end of order 2


  % fill the symmetric components of g
  for i=1:4
    for j=i+1:4
      g(:,:,:,j,i) = g(:,:,:,i,j);
    end
  end
  





    % all of the functions expect a single number for T, and a vector of z at a fixed point (x,y)

    % *********************
    %   h - z^2
    % *********************

    function resh2s6 = h2s6(z, T) % 1/18/13
      % has z.^2 factored out. but only z, no T (that goes for everybody)!!
      % -1/2r^4
      x = z(:).*T;
      resh2s6 = 0.*z(:) +  (T.^2)./4;
    end

    function resh2s7 = h2s7(z, T)
      % has z.^2 factored out. but only z, no T!!
      %resh2b = 0;
      x = z(:).*T;
      if params.metricLinearFnInterpolation
        resh2s7 = myinterp1(interpgrid, nh2s7, x);
      else
        resh2s7 = cheb_interp(nh2s7, x, 0, 2);
      end
      resh2s7 = resh2s7 .* T .* T;
    end

    % *********************
    %   k - 1/z^2
    % *********************

    function resk2s3 = k2s3(z, T) % 1/18/13
      % k2s3 with 1/z^2 extracted
      % 2 S3
      resk2s3 = 0.*z(:) - 1;
      resk2s3 = resk2s3./(T.^2);
    end
    function resk2s4 = k2s4(z, T)
      % placeholder
      resk2s4 = 0.*z(:);
      resk2s4 = resk2s4./(T.^2);
    end
    function resk2s5 = k2s5(z, T) % 1/18/13
      % k2s5 with 1/z^2 extracted
      % 1/2 S5
      resk2s5 = 0.*z(:) - 1./4;
      resk2s5 = resk2s5./(T.^2);
    end
    function resk2s6 = k2s6(z, T) % 1/18/13
      % k2s6 with 1/z^2 extracted.
      % -(1+4r^3)/(2r^3) S6
      x = z(:).*T;
      resk2s6 = 1+(-1./8).*x.^3;
      resk2s6 = resk2s6./(T.^2);
    end
    function resk2s7 = k2s7(z, T)
      % k2s7 with 1/z^2 extracted.
      x = z(:).*T;
      if params.metricLinearFnInterpolation
        resk2s7 = myinterp1(interpgrid, nk2s7, x);
      else
        resk2s7 = cheb_interp(nk2s7, x, 0, 2);
      end
      resk2s7 = resk2s7./(T.^2);
    end

    % *********************
    %   j - 1/z
    % *********************

    function resj2v3 = j2v3(z, T) % 1/18/13
      % j2v3 with 1/z extracted!
      % same as in mark's paper
      % 1/(2r^2) V3
      resj2v3 = 0.*z(:) + -1./2;
      resj2v3 = resj2v3./T;
    end

    function resj2v4 = j2v4(z, T) % 1/18/13, same in march 03/02/13
      % with 1/z extracted
      % same as in mark's paper!
      % -1/(2r^2(1+r+r^2)) V4
      resj2v4 = zeros(length(z), 1);
      for i=1:length(z)
        x = z(i).*T;
        if (x < 0.001)
          resj2v4(i) = (-1./4).*x.^2+(1./10).*x.^3+(-1./28).*x.^5+(1./40).*x.^6+(-1./70).*x.^8+(1./88).*x.^9+(-1./130).*x.^11+(1./154).*x.^12+(-1./208).*x.^14+(1./238).*x.^15+(-1./304).*x.^17+(1./340).*x.^18+(-1./418).*x.^20; % feb 27
        else
          resj2v4(i) = (1./54).*x.^(-2).*(9.*((-2)+x).*x+2.*3.^(1./2).*pi.*((-1)+x.^3)+(-12).*3.^(1./2).*((-1)+x.^3).*atan(3.^(-1./2).*(1+2.*x)));
        end
      end
      resj2v4 = resj2v4./T;
    end

    function resj2v5 = j2v5(z, T) % 1/18/13
      %j2v5 with 1/z extracted, the resulting thing has lowest order const.
      % (-1+r+3r^2+2r^3+r^4)/(2r^2(1+r+r^2)^2) (V5-V6)
      resj2v5 = zeros(length(z), 1);
      for i=1:length(z)
        x = z(i).*T;
        if (x < 0.001)
          resj2v5(i) = (-1./2)+(-1./2).*x.^2+(1./5).*x.^3+(-1./14).*x.^5+(1./20).*x.^6+(-1./35).*x.^8+(1./44).*x.^9+(-1./65).*x.^11+(1./77).*x.^12+(-1./104).*x.^14+(1./119).*x.^15+(-1./152).*x.^17+(1./170).*x.^18+(-1./209).*x.^20;
        else
          resj2v5(i) = real((-1./54).*x.^(-2).*(36.*x+9.*x.^2+4.*3.^(1./2).*(pi+(-6).*atan(3.^(-1./2).*(1+2.*x)))+4.*3.^(1./2).*x.^3.*(2.*pi+(sqrt(-1).*3).*log(sqrt(-1)+(-1).*3.^(1./2)+(sqrt(-1).*2).*x)+(sqrt(-1).*(-3)).*log(sqrt(-1)+3.^(1./2)+(sqrt(-1).*2).*x))));
        end
      end
      resj2v5 = resj2v5./T;
    end

    function resj2v6 = j2v6(z, T) 
      resj2v6 = -1.*j2v5(z,T);
    end

    function resj2v7 = j2v7(z, T) % 3/02/13
      resj2v7 = z.*0;
    end




    % ***** % 
    %   a   %
    % ***** %

    % no extraction

    function resa2t2 = a2t2(z, T)
      x = z.*T;
      if params.metricLinearFnInterpolation
        resa2t2 = myinterp1(interpgrid, na2t2, x);
      else
        resa2t2 = cheb_interp(na2t2, x, 0, 2);
      end
    end

    function resa2t3 = a2t3(z, T)
      resa2t3 = a2t2(z,T);
    end

    function resa2t4 = a2t4(z, T)
      x = z.*T;
      if params.metricLinearFnInterpolation
        resa2t4 = myinterp1(interpgrid, na2t4, x);
      else
        resa2t4 = cheb_interp(na2t4, x, 0, 2);
      end
    end
    function resa2t5 = a2t5(z, T)
      x = z.*T;
      if params.metricLinearFnInterpolation
        resa2t5 = myinterp1(interpgrid, na2t5, x);
      else
        resa2t5 = cheb_interp(na2t5, x, 0, 2);
      end
    end


end
